##############################
# STATE REACH AND DEVELOPMENT // PSRM // Replication
#
# Voronoi cells vs. square grid cells comparison
#
##############################


# GLOBALS #######################

# Cell sizes
vor.cell.size <- c(100, 400, 1600, 6400)


# FUNCTIONS #######################
source("code/voronoi_functions.R")


# DATA ############################

# Locals: Data for Nigeria
cowcode <- 475
year <- 1992:2013

# Load Uganda Regions ##
admin.shp <- readOGR(dsn = file.path("data", "africa_regions_panel.GeoJSON"), 
                     layer="africa_regions_panel")
admin.shp <- admin.shp[admin.shp$startyr <= max(year) & admin.shp$endyr >= min(year) & admin.shp$cowcode == cowcode,]

# Compute all constitutive part
adm.const.parts <- get_adm_const_parts(admin.shp = admin.shp, base.year = 1992, ncore = 10, size.cutoff = .001)
adm.const.parts$poly.id <- 1:length(adm.const.parts)
adm.const.parts <- gBuffer(adm.const.parts, byid = T, width = 0) 
adm.const.parts$adm.id <- 1:length(adm.const.parts)

# MAKE CELLS ##################

# Voronoi Cells ###

# Draw voronoi units
voronoi.ls <- lapply(vor.cell.size, function(i){
  print(paste("Make voronoi cells of size", i, "km2"))
  sample_vorcells(spdf = adm.const.parts, res = .01, size = i, size.km2 = T, 
                  seed = 1, iter.max = 400,
                  sample.type = "nonaligned", ncore = 10)
})

# Add cowcodes
voronoi.ls <- lapply(voronoi.ls, function(v){
  v$adm.id <- adm.const.parts$adm.id[v$unit.id]
  v
})

# Grid Cells ###

# Master Grids
grids.ls <- lapply(vor.cell.size, function(s){
  ext <- extent(admin.shp)
  res <- ((s^.5)/111.32)
  cols <- ceiling((ext[2] - ext[1])/res)
  rows <- ceiling((ext[4] - ext[3])/res)
  raster(xmn=ext[1], xmx=ext[1] + cols*res, ymn=ext[3], ymx=ext[3] + rows*res, 
         resolution = res, vals = NA)
})

# Rasters to polygon cells
cell.ls <- lapply(grids.ls, function(r){
  rasterToPolygons(r, na.rm = F)
})

# Crop by countries

# ... setup cluster
cl <- makeCluster(getOption("cl.cores", ncore))
clusterExport(cl, list("cell.ls", "adm.const.parts"))
registerDoParallel(cl)

# ... crop cells
cell.crop.ls <- lapply(1:length(cell.ls), function(i){
  clusterExport(cl, list("i"), envir = environment())
  foreach(c = 1:length(adm.const.parts),
          .packages = c("raster"), .noexport = c("cell.ls", "adm.const.parts"),
          .options.multicore = list(preschedule = FALSE), .combine = rbind) %dopar% {
            r <- raster::crop(cell.ls[[i]], adm.const.parts[c,])
            r$adm.id <- adm.const.parts$adm.id[c]
            r
          }
})

# ... stop cluster
stopCluster(cl)


# EVALUATION ##################

# FIGURE A10: Simple comparison plot

# ... prepare
plot.size <- which(vor.cell.size == 400)

# ... plot Grid
png(file.path(fig.path, paste0("figa10_vor_vs_cell_ex1", ".png")), 
    height = 4, width = 4, res = 300, unit = "in")
par(mar = c(0,0,2,0))
plot(cell.crop.ls[[plot.size]],
     border = "black", main = "Nigeria, regular grid, \n 400 sq. km", cex = .5)
dev.off()

# ... plot voronoi
png(file.path(fig.path, paste0("figa10_vor_vs_cell_ex2", ".png")), 
    height = 4, width = 4, res = 300, unit = "in")
par(mar = c(0,0,2,0))
plot(voronoi.ls[[plot.size]],
     border = "black", main = "Nigeria, Voronoi cells, \n 400 sq. km", cex = .5)
dev.off()


# FIGURE A11: Size distributions

# ... mean size admin units
adm.size <- areaPolygon(adm.const.parts)/1e6
mean(adm.size)
median(adm.size)

# ... Prepare data
plot.df <- rbind(data.frame(type = "Voronoi",
                            size = unlist(lapply(voronoi.ls, areaPolygon)),
                            area = rep(vor.cell.size, unlist(lapply(voronoi.ls, length)))),
                 data.frame(type = "Regular Grid",
                            size = unlist(lapply(cell.crop.ls, areaPolygon)),
                            area = rep(vor.cell.size, unlist(lapply(cell.crop.ls, length)))))
plot.df$size <- plot.df$size / 1e6
plot.df$target <- plot.df$area
plot.df$area <- paste0(plot.df$area, " square km")
plot.df$area <- factor(plot.df$area, levels = unique(plot.df$area))

# ... sumstats
sumstat.labs <- rbind(aggregate.data.frame(list(value = plot.df$size),
                                           plot.df[,c("area", "type")],
                                           FUN = sd),
                      aggregate.data.frame(list(value = plot.df$size),
                                           plot.df[,c("area", "type")],
                                           FUN = sd),
                      aggregate.data.frame(list(value = plot.df$size),
                                           plot.df[,c("area", "type")],
                                           FUN = mean),
                      aggregate.data.frame(list(value = plot.df$size),
                                           plot.df[,c("area", "type")],
                                           FUN = median))
sumstat.labs$value <- as.character(signif(sumstat.labs$value , 3))
sumstat.labs$stat <- rep(c("Title","SD", "Mean", "Median"), each = length(vor.cell.size) * 2)
sumstat.labs$value[sumstat.labs$stat == "Title"] <- ""
sumstat.labs$stat[sumstat.labs$stat == "Title"] <- rep(c("Voronoi", "Reg. Grid"), each = length(vor.cell.size))
sumstat.labs$label <- paste0(sumstat.labs$stat, ":", " ", sumstat.labs$value)

# ... sumstats positions
max.x <- unlist(lapply(unique(plot.df$target), function(s){
  max(unlist(lapply(unique(plot.df$type), function(t){
    density(plot.df$size[plot.df$type == t & plot.df$target == s])$y
  })))
}))
sumstat.labs$pos.y <-  unlist(lapply(1:4, function(y){
  rep(unlist(lapply(max.x, function(x){seq(x, 0, length.out = 15)[y]})), 2)}))
sumstat.labs$pos.y <- sumstat.labs$pos.y  -  c(rep(0, length(vor.cell.size)), max.x / 3)
sumstat.labs$pos.x <- rep(c(0), each = 4)

# ... plot
png(file.path(fig.path, paste0("figa11_vor_vs_cell_size", ".png")), 
    height = 4, width = 10, res = 300, unit = "in")
print(ggplot(plot.df, aes(x = size, col = type, lty = type, group = type)) +
  geom_line(stat = "density", lwd = 1) + 
  xlab(~"Size of units " (km^2)) +
  geom_vline(data = unique(plot.df[,c("area", "target")]), aes(xintercept = target), 
             lty = 2, col = "darkgrey") +
  theme_minimal() +
  geom_text(data = sumstat.labs, aes(x = pos.x, y = pos.y, label = label),
            hjust = 0, cex = 2, col = "black") +
  facet_wrap(~ area, nrow = 1, scale = "free")) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 3)) 
dev.off()


# FIGURE A14: Compactness: Polsby-Popper Test ###
# https://en.wikipedia.org/wiki/Polsby-Popper_Test

# Function
calc_ppt <- function(spdf){
  area <- areaPolygon(spdf)/1e6
  perimeter <- perimeter(spdf)/1e3
  (4*pi * area)/(perimeter^2)
}

# Voronoi
vor.ppt <- lapply(voronoi.ls, calc_ppt)

# Grid
cell.ppt <- lapply(cell.crop.ls, calc_ppt)


# Plot

# ... Prepare data
plot.df <- rbind(data.frame(type = "Voronoi",
                            pps = unlist(vor.ppt),
                            area = rep(vor.cell.size, unlist(lapply(vor.ppt, length)))),
                 data.frame(type = "Regular Grid",
                            pps = unlist(cell.ppt),
                            area = rep(vor.cell.size, unlist(lapply(cell.ppt, length)))))
plot.df$area <- paste0(plot.df$area, " square km")
plot.df$area <- factor(plot.df$area, levels = unique(plot.df$area))

# ... sumstats
sumstat.labs <- rbind(aggregate.data.frame(list(value = plot.df$pps),
                                           plot.df[,c("area", "type")],
                                           FUN = sd),
                      aggregate.data.frame(list(value = plot.df$pps),
                                           plot.df[,c("area", "type")],
                                           FUN = sd),
                      aggregate.data.frame(list(value = plot.df$pps),
                                          plot.df[,c("area", "type")],
                                          FUN = mean),
                     aggregate.data.frame(list(value = plot.df$pps),
                                          plot.df[,c("area", "type")],
                                          FUN = median))
sumstat.labs$value <- as.character(signif(sumstat.labs$value , 3))
sumstat.labs$stat <- rep(c("Title","SD", "Mean", "Median"), each = 2 * length(vor.cell.size))
sumstat.labs$value[sumstat.labs$stat == "Title"] <- ""
sumstat.labs$stat[sumstat.labs$stat == "Title"] <- rep(c("Voronoi", "Reg. Grid"), each = length(vor.cell.size))
sumstat.labs$label <- paste0(sumstat.labs$stat, ":", " ", sumstat.labs$value)

# ... sumstats positions
sumstat.labs$pos.y <- rep(25 - c(1:4)*2, each = 2 * length(vor.cell.size)) - 
  rep(c(0, 10), each = length(vor.cell.size))
sumstat.labs$pos.x <- rep(c(0), each = 4)

# ... plot
png(file.path(fig.path, paste0("figa12_vor_vs_cell_pps", ".png")), 
    height = 4, width = 10, res = 300, unit = "in")
print(ggplot(plot.df, aes(x = pps, col = type, lty = type, group = type)) +
  geom_line(stat = "density", lwd = .8) + 
  xlab("Polsby-Popper (compactness) score of units") +
  coord_cartesian(ylim=c(0, 25), xlim=c(0,1)) +
  theme_minimal() +
  geom_text(data = sumstat.labs, aes(x = pos.x, y = pos.y, label = label),
            hjust = 0, cex = 2, col = "black") + 
  scale_x_continuous(breaks = seq(0,1, by = .5), 
                     labels = c(0, .5, 1)) +
  facet_wrap(~ area, nrow = 1))
dev.off()



# Clean up
rm(list = c("adm.const.parts", "admin.shp", "cell.crop.ls", 
            "cell.ls", "grids.ls", "voronoi.ls", "cl", "plot.df"))


